function [output] = IM_main(input_array)
% This function uses a sequence of clustering structure to perform 
% Ibragimov and Muller (2010). It produces a vector of test restuls of 0 or
% 1. "testIM_star" is the test result using the clustering selected 
% according to a size-power tradeoff criterion. 

% data input
method = input_array.method;
data = input_array.data;
Sigma_hat = input_array.Sigma_hat;
group_matrix = input_array.group_matrix_km;
alpha_sig = input_array.alpha_sig;
alt_vec = input_array.alt_vec;
alt_vec_tradeoff = input_array.alt_vec_tradeoff;

G_vec = max(group_matrix); % vector of group numbers
l_g = length(G_vec);

% t1_vec = zeros(1,length(G_vec));
ave_power_vec = zeros(1,length(G_vec));
p_adj_vec = zeros(1,length(G_vec));
rej_rate_vec = zeros(1,length(G_vec));
% adjust p-value thresholds to control sizes for each clustering
for jj = 1 : length(G_vec)
    membership = group_matrix(:,jj);

    [rej_rate,p_adj] = rej_rate_IM_mean(Sigma_hat,membership,alpha_sig,0);
    rej_rate_vec(jj) = min(rej_rate,alpha_sig);
    p_adj = min(p_adj,alpha_sig); % cannot be more aggressive

%     type_I = type1IM_mean(Sigma_hat,membership,p_adj);
%     ave_power_vec(jj) = type_I;

    power_vec = zeros(length(alt_vec_tradeoff),1);
    for i_alt = 1 : length(alt_vec_tradeoff)
        power_vec(i_alt) =  rej_rate_IM_mean(Sigma_hat,...
            membership,p_adj,alt_vec_tradeoff(i_alt));
    end
    ave_power = mean(power_vec);
    
%     lossVec(jj) = type_II;
%     t1_vec(jj) = type_I;
    ave_power_vec(jj) = ave_power;
    p_adj_vec(jj) = p_adj;
end
[ave_power,ind_G] = max(ave_power_vec);
group_star = group_matrix(:,ind_G);
p_adj = p_adj_vec(ind_G);
G_star = G_vec(ind_G);
estimated_size = rej_rate_vec(ind_G);


%% ESTIMATION
D = data.D;
X = data.X;
Y = data.Y;
        
switch lower(method)
    case 'ols'
        Z = D;
        
        alt_power_curve = linspace(alt_vec(1),alt_vec(end),50);

    case 'iv'
        Z = data.Z;
        
        alt_power_curve = linspace(alt_vec(2),alt_vec(end-1),50);
        
end

[bfmtmp,sfmtmp] = FamaMacbeth(D,X,Y,Z,group_star); % FM results
thetaHat_fm = bfmtmp(1);
se_fm = sfmtmp(1);


%% OUTPUT for estimation and inference table
% p-value under the null
t_stat = (thetaHat_fm-0)/se_fm;
p_null = 2*(1-tcdf(abs(t_stat),max(group_star)-1));

% p-value under the alternatives
t_stat_alt = (thetaHat_fm+alt_vec-0)/se_fm;
p_alt = 2*(1-tcdf(abs(t_stat_alt),max(group_star)-1));

output.estimate = thetaHat_fm;
output.p_adj = p_adj;
output.p_null = p_null;
output.p_alt = p_alt;


%% OUTPUT for power curve
% alt_power_curve = linspace(alt_vec(1),alt_vec(end),50);
t_stat_power_curve = (thetaHat_fm+alt_power_curve-0)/se_fm;
p_power_curve = 2*(1-tcdf(abs(t_stat_power_curve),max(group_star)-1));

output.alt_power_curve = alt_power_curve;
output.p_power_curve = p_power_curve;


%% p-value with fixed G

p_fixed_G = zeros(1,length(G_vec));

for i_g = 1 : l_g
    
group = group_matrix(:,i_g);

[bfmtmp,sfmtmp] = FamaMacbeth(D,X,Y,Z,group); % FM results
thetaHat_fm_temp = bfmtmp(1);
se_fm = sfmtmp(1);

% p-value under the null
t_stat = (thetaHat_fm_temp-0)/se_fm;
p_fixed_G(i_g) = 2*(1-tcdf(abs(t_stat),max(group)-1));

end


%% OUTPUT for clustering table
output.G_star = G_star;
output.ave_power = ave_power;
output.p_fixed_G = p_fixed_G;
output.p_threshold_fixed_G = p_adj_vec;
output.estimated_size = estimated_size;
% p-value used for comparing actual power and estimated power
t_stat_tradeoff = (thetaHat_fm+alt_vec_tradeoff-0)/se_fm;
p_tradeoff = 2*(1-tcdf(abs(t_stat_tradeoff),max(group_star)-1));
output.p_tradeoff = p_tradeoff;

